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Abstract 

We evaluate by Monte Carlo simulations various singular thermodynamic 
quantities X, for an ensemble of quenched random Ising and Ashkin- Teller 
models. The measurements are taken at T c and we study how the distribu- 
tions P(X) (and, in particular, their relative squared width, Rx) over the 
ensemble depend on the system size I. The Ashkin- Teller model was stud- 
ied in the regime where bond randomness is irrelevant and we found weak 
self averaging; Rx ~ l~ — > 0, where a < and v are the exponents (of the 
pure model fixed point) governing the transition. For the site dilute Ising 
model on a cubic lattice, known to be governed by a random fixed point, 
we find that Rx tends to a universal constant independent of the amount of 
dilution (no self averaging). However this constant is different for canonical 
and grand canonical disorder. We identify the pseudo-critical temperatures 
T c (i, I) of each sample i, as the temperature at which the susceptibility reaches 
its maximal value. The distribution of these T c (i,l) over the ensemble was 
investigated; we found that its variance scales as (5T c (l)) 2 ~ l~~ . These 
results are in agreement with the recent predictions of Aharony and Harris. 
Our previously proposed finite size scaling ansatz for disordered systems was 
tested and found to hold. When we fit the data obtained for many samples 
of different sizes by a sample-independent form, the resulting scaling function 
was found to be universal and to behave similarly to pure systems. We did 
observe that to describe deviations from this universal function we do need 
sample-dependent scaling functions. These deviations are, however, relatively 
small and this led us to an interesting side result: sample-to-sample fluctua- 
tions of x max , the susceptibility measured at T c (i, I), are smaller by a factor of 
70 than those of x(T c )- This indicates that to obtain a fixed statistical error 
it may be more computationally efficient to measure ^ max . 
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I. INTRODUCTION 



How is the critical behavior affected by the introduction of disorder (usually dilution 
or bond-randomness) into a model? This question has been extensively studied [|l| experi- 
mentally, analytically |§ and numerically || for quite some time now. The Harris criterion 
IIHH states that 0, the scaling index of the operator corresponding to randomness at the 
pure system fixed point (also called the crossover exponent ) is equal to 2e (« p and v p are 
the specific- heat and correlation length exponents of the pure model). Thus the critical 
behavior of the pure system (p) is unaltered by disorder if a p < and a weakly disordered 
system will have the same critical exponents as the pure one. If a p > even a weakly 
disordered system will not belong to the same universality class as the pure one. If a p — 
the situation is marginal. 

In Renormalization Group (RG) calculations it is possible to determine to which fixed 
point a certain disordered model flows and to determine the nature of this fixed point. A 
disordered fixed point is characterized by a fixed distribution (of finite width) of couplings 
while a pure fixed point is characterized by a delta function type distribution: a single 
coupling set. On the other hand, to the best of our knowledge, in the current Monte Carlo 
state of the art, there is no method that can check directly whether a certain model is 
governed by a pure or disordered fixed point. Most Monte Carlo studies concentrated on 
calculating critical exponents of a certain disordered model. If these exponents agreed with 
an RG calculation, this served as indirect evidence to the nature of the governing fixed 
point. In such numerical and experimental studies of disordered systems near their critical 
point finite samples with quenched disorder are used; any sample i is a particular random 
realization of the quenched disorder. A measurement (or calculation) of any density of an 
extensive thermodynamic property X (e.g. X = E, M, Ch or x) would yield a different value 
Xi for every sample i because of the differences in the realizations of the quenched disorder. 
Here Xi is the exact thermal average of the sample, which in an experiment or a numerical 
study can only be estimated by X*. In an ensemble of disordered samples of linear size / the 
values of Xi are governed by a probability distribution P(X). In most MC studies only the 
ensemble average [X] is studied. As is shown in this study, it is possible to obtain direct 
evidence, by MC, as to the nature of the governing fixed point by studying P(X) and the 
factors which govern its shape. As it turns out, this can be done by studying the question 
of self averaging, which concerns the behavior of the width of P{X) as the system size is 
increased. We characterize P(X) by the ensemble average [X] and variance 

V X = [X 2 ] - [X? . (1) 

Suppose that X is a singular density of an extensive thermodynamic property, such as 
M, x or the singular part of E or C . The system is said to exhibit self- averaging || if 

R x (l) = V X /[X} 2 -4fl as Noo, (2) 

otherwise, if Rx tends to a constant different from zero, it is said to exhibit lack of self 
averaging. In a self averaging system a single very large sample is a sufficient representative 
of the ensemble. But in a non self-averaging system a measurement performed on a single 
sample, no matter how large, does not give a meaningful result and therefore must be 
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repeated using many samples from the ensemble. In a MC study of a self averaging disordered 
system the number of samples needed to obtain [X] to a given relative accuracy decreases 
with increasing I. On the other hand, in a non self averaging system the number of samples 
which must be simulated is independent of I and the total amount of work rises very strongly 
with I. 

Off criticality, where I is much larger than the correlation length £, as first argued by 
Brout J7|, we may divide the sample % into n large subsamples (much larger than £). If 
we assume that the coupling between neighboring subsystems is negligible, then the value 
of any density of an extensive quantity over the whole sample X^ is equal to the average 
of the (independent) values of this quantity over the subsamples. Provided the probability 
distribution of the X's of the subsamples has a finite variance, then according to the Central 
Limit theorem the value of Xi is distributed with a Gaussian probability distribution around 
its mean [Xi\. The square of the width of the Gaussian, Vx, is proportional to - ~ l~ d . In 
this case (§) is fulfilled, and X is called self-averaging. In such a case, where Rx ~ l~ d , X 
is called strongly self averaging ||. 

Close to criticality, where £ ~ I, the Brout argument does not hold, since the average of 
X over neighboring subsamples may not be considered as independent. Thus at criticality 
there is no reason to expect that Rx ~ l~ d - In a previous study we considered the 
question of self-averaging at the critical temperature of the infinite lattice, T c °°, through 
Monte Carlo simulations of several random-bond Ashkin Teller models. Our findings, which 
are summarized shortly below, prompted us to suggest a heuristic finite size scaling theory 
for disordered systems, based on physical considerations similar to those leading to the 
Harris criterion. We characterized every sample i of size I by a pseudo-critical temperature 
T c (i, I) and introduced the reduced temperature of each sample i 

T — TJi, I) , , 

U = T^ 1 ■ (3) 

J- c 

T c (i,l) was assumed to fluctuate around its ensemble average T c (l) with width 5T c (l). We 
then assumed [for samples with T close to T c (i, I)] a sample dependent finite size scaling 
form 

X i (T,l)=FQ i (t i l m ). (4) 

Here p is the exponent characterizing the behavior of [X] at T c ; e.g. for X = x, P = ^ • The 
form of the function Qi(Z) (or its coefficients) was assumed to be sample dependent but 
the critical exponents p and yt — - were assumed to be universal or sample independent. 
Relying on (|) we then related Rx at T c °° to 8T c (l). By assuming that 

(6T c (l)f~l- d (5) 

the theory predicted that when the specific heat exponent of the disordered model a < 
then, at criticality, 

R x ~l» . (6) 
Thus, if - = 0, X is non self-averaging, but if — d < - < 0, X is called weakly self averaging 



Note that according to [10| - < in any disordered system, though claims to the 



contrary exist (e.g. |TTJ and [pqj ). 
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In a subsequent study, Aharony and Harris |13j (AH) used renormalization group analysis 
to study the dependence of P{X) on I and £. For I ^> £ they recovered strong self-averaging: 
P(X) approaches a Gaussian with relative variance 

Rx ~ (l/0' d - For Z < £ they found 
two different behaviors. When randomness was irrelevant and the system was governed by 
a pure fixed point they found 

R x ~ . (7) 

In this case the critical exponents of the disordered model are the same as those of a pure 
model, ~ = (2^ ) so that (0) is in agreement with @. On the other hand, when the 

system is governed by a random fixed point, they found that P{X) approaches a universal, I 
independent shape, and Rx — > const as / — > oo, which implies lack of self-averaging. When 
a of the random model is negative, this prediction is in disagreement with (|6]). As AH point 
out, this disagreement between the RG result and our scaling ansatz can be reconciled if we 
assume that in disordered systems governed by a random fixed point, (||) is substituted by 

{5T c {l)f - I' 1 ■ (8) 

In the Monte Carlo study of || several bond-disordered Ashkin- Teller models on a square 
lattice were simulated. These included the bond-disordered Ising model where (jf) = 

+ (C ~ log/ at T c °°), and several bond-disordered Ashkin- Teller models where f^J 



> 



0. According to renormalization group expansions |l^,p|JT5[| around the pure Ising model 
these models are governed by the pure Ising fixed point, while according to MC results 
fl6|| the bond-disordered Ashkin- Teller models may be governed by a different (possibly 
random) fixed point with - = + . It was found that far from criticality all thermodynamic 
quantities which were examined (energy, magnetization, specific heat, susceptibility) are 
strongly self averaging, that is Rx ~ l~ d - At criticality though, the results indicated that 
the magnetization m and the susceptibility \ are non self averaging. The energy E at 
criticality was clearly weakly self averaging, that is Ve ~ l~ Vv with < y v < d (Here E 
includes the analytic and singular parts of the energy). The theory (^) is not applicable in 
the asymptotic limit (I — > oo) to the bond-disordered Ashkin- Teller model where ^ = + . 
Nonetheless in the accessible range of lattice sizes good agreement between the theory and 
the data for V x and Ve was found. In particular R x was shown to behave very similarly to 
the specific heat C, as suggested by (||), for a wide range of lattice sizes and for different 
degrees of disorder. In a very recent MC study of a mean field Potts glass a lack of self 



averaging of the order parameter was found as well |17 . 

In this paper we set out to resolve three issues, neither of which could be investigated 
by studying the Ashkin Teller model at a p > 0. 

1. So far the prediction (0) has not been tested numerically for the case a p < 0. In this 
case randomness is irrelevant, a = a p and the AH prediction (|7]) coincides with (|6]). 

2. When a p > but at the random fixed point a < 0, the scaling theory predicts 
weak self averaging, in disagreement with the AH result Rx — > const. Resolution of 
this discrepancy will shed light on whether the ansatz (0) or the AH prediction ([Bp 
governs the width of the distribution of the pseudocritical temperatures ST c (l). 
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3. Testing the scaling form (|]) : We wish to determine whether it holds and whether the 
sample dependence enters only via ij or also through the scaling function Qi. 

With these three goals in mind we set out to examine by Monte Carlo simulations the 
question of self averaging at criticality in two different models. The first is a bond-disordered 
Ashkin- Teller model at a point where (j^j « —0.54 (a large negative value was chosen to 

yield unambiguous results). To address the second, more important issue, we simulated 
the site-dilute Ising model on a cubic lattice. Because the critical exponent a p is positive 



for the three dimensional pure Ising model, a p ~ 0.11 | IBfl , the Harris criterion predicts 
that randomness will lead to a different critical behavior. The disordered model is believed 
19| , |20[| to be governed by a random fixed point with a < 0. Thus, according to AH, this 



model should exhibit lack of self averaging while according to our finite size scaling theory 
it should exhibit weak self averaging as in (P), if the assumption @ is valid. 

Besides calculating Rx at the critical point, we decided to test directly which one of @ 
or (H) is correct. To this end we calculated the pseudo-critical temperature of each sample 
taken from an ensemble of site-diluted 3D Ising samples at different lattice sizes. The pseudo- 
critical temperature was defined as the temperature of the maxima of the susceptibility of 
that sample. This calculation was done by using the histogram reweighting method |2T|-|2S[. 
This method allows to use the results of a simulation at one temperature for calculating 
the value of thermodynamic observables in a whole range of nearby temperatures. We thus 
obtained for each lattice size distributions of pseudo-critical temperatures T c (i,l) and were 
able to calculate their mean, T c (l), and variance, (5T c (l)) 2 . 

To investigate the extent to which the scaling form (|4]) holds, we studied the relationship 
between the sample dependent magnetization mj(T c ) and T c (i,l) using the data collapse 
technique. We did find convincing support for the finite size scaling ansatz (|^) but also 
found evidence for sample dependence of the scaling function . 

This work is organized as follows. In the first part of Sec. [TI] we define the random bond 
Ashkin- Teller model which was simulated and summarize its critical behavior as found by 
the finite size scaling results. In the second part of Sec. [TT] we give our results concerning 
self-averaging at criticality. The results indicate clearly that Rx is weakly self averaging and 



are in good agreement with (|7|). In Sec. ( I I] we summarize some relevant properties of the 
site dilute Ising model on a cubic lattice and give some details of the simulation. Finite size 
scaling results for some observables at criticality are given as well. In section [V] we analyze 
and discuss our results concerning self averaging at T£P. These results seem to indicate the 
correctness of the AH scenario, whereby Rx is non self averaging. In Sec. [V| we study the 
distributions of the pseudo-critical temperatures of the site dilute Ising model. The scaling 
of (<5T C (Z)) 2 does not agree with @ but seems to agree with (|8|), giving additional evidence 
for the validity of the AH scenario. In Sec. [V| we also analyze the distributions of the 
maximal susceptibilities x{T c {h 0)) an d investigate the extent to which the scaling form (Eh 



holds. The work is summarized in Sec. VI 



II. WEAK SELF AVERAGING IN AN ASHKIN-TELLER MODEL WITH 

IRRELEVANT DISORDER 
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A. Definition of the model 



The model we study is the random-bond Ashkin- Teller model on a square lattice. On 
every site of the lattice two Ising spin variables, <7j and Tj, are placed. Denoting by (ij) a 
pair of nearest neighbor sites, the Hamiltonian is given by 

H = - Y^[Ki,j{(Ti(Tj + TiTj) + AijCTiTiajTj] . (9) 

(ij) 

K it j and Ay are chosen according to 

, x _ / (A' 1 , A 1 ) with probability \ , , 

I a m; a mJ - | (K 2 ,A 2 ) with probability § " 1 UJ 

The homogeneous (or pure) model [^6| [(A 1 , A 1 ) = (if 2 , A 2 ) = (if , A)] possesses a line 
of critical points along which critical exponents vary continuously, so that it flows under 
Renormalization Group (RG) onto a line of fixed points. The scaling exponent corresponding 
to randomness, <p = (a/v), which is analytically known fl27U28H , also varies continuously 
along this line. However along the part of this line ( A > 0) interpolating between the Ising 
(A = 0) and four state Potts (A = if) models it takes positive values, (1 > > 0), so that 
randomness is relevant. Indeed the critical behavior of the disordered model jT5|JT^ ,|^JT6|] 
was found to be different from that of the pure one. The self averaging properties of the 
disordered model were examined in ||, and as discussed in the introduction, a lack of self 
averaging was found. Along the other part of this line ( A < 0) = ^ is negative, so that 
according to the Harris criterion randomness is irrelevant. Thus a slightly disordered model 
will exhibit the same behavior as its pure version. It turns out, according to this study 
(in agreement with RG calculations p,15[|) that a model with finite disorder will flow under 



RG onto a pure model that is along the part of the line of fixed points where </>=-< 0. 
Thus according to AH as well as according to (|Bp we expect to find weak self-averaging at 
criticality for disordered versions of models corresponding to this part of the Ashkin Teller 
critical line (A < 0). 

The Ashkin Teller model is a convenient paradigm for studying critical behavior of disor- 
dered systems both because its scaling exponent corresponding to randomness <p — ^ varies 
continuously and because part of the critical manifold of the random-bond model can be 
found exactly through duality. In parts of the coupling space where only two phases exist 
the self dual manifold 

(if 2 ,A 2 ) = (if\Ai) (11) 

is critical. Here (if 1 , A 1 ) are the dual couplings of (if 1 , A 1 ) according to the duality trans- 
formation of the Ashkin- Teller model; a discussion of this point can be found in ||16|| . Since 



the extent of deviation from pure behavior is obviously determined by the difference between 
the two sets of couplings, we have chosen to study a model with the ratio 

if 2 1 

T^To- < 12 > 

so that randomness will be pronounced. In addition, a ratio of 
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A 1 9 

iP = -f = -v> (13) 

was chosen. Since - decreases as / increases and - = for / = ( Ising model), we have 
chosen / = yjj so that ^ will be a pronounced negative number. Equations (|lll - |l3| ) define the 
couplings of the model simulated, where the temperature T was absorbed into the couplings 
Kij,Aij. 

The critical behavior of the disordered model is compared with that of an Anisotropic 
Ashkin-Teller model which is used as a reference pure model [[TJJ. This model has the same 
Hamiltonian (^) but with the couplings distributed as follows: 

(K A ) — < ^ or bonds (hj) m the horizontal direction , . 

1,3 ) (K 2 ,A 2 ) for bonds in the vertical direction 



In the Monte Carlo simulations we used a cluster algorithm [^] which is described in 



TM. The main idea is to embed into the Ashkin-Teller model an Ising model and simulate 



it using the Wolff |30| single cluster algorithm for the Ising model. The number of samples 
simulated was n = 2000 for linear lattice sizes 4 < I < 64, n(128) = 1200 for I = 128, 
and n(256) = 436. For each sample i Monte Carlo estimates of various observables Xi and 
their errors 5Xi were calculated. Next we list results concerning the critical behavior of the 
estimated bond-disordered ensemble averages [Xj\ as a function of lattice size. 



B. Critical behavior of the model 

Here we give the finite size dependences of averages (over all samples) of various observ- 



ables, defined as in [16], at the critical point T c °° defined through ( 1 1 -p"3|) . 

In figure [l] we plot the specific heat of the random bond and the anisotropic models as 
a function of log/. The solid lines are fits to the finite size scaling form 

C = B Q + B x lv . (15) 

Using lattice sizes of 16 < / < 256, we find - = —.745(4) for the anisotropic model, while 
using lattice sizes of 24 < I < 256 , we find - = —.536(32) for the random bond model 
(note: B\ is negative). Thus this strongly disordered model apparently flows under RG onto 
a pure model with different exponents than its anisotropic version but still one that is along 
the part of the line of fixed points where < 0. 

For both models the magnetization m and susceptibility \ were fitted to the forms 

m = A m l v ^ = A x l» . (16) 

Similarly the polarization P (magnetization of the r spins) and susceptibility of the polar- 
ization x were fitted to the forms 

P = Ap l-— x (p) = A xW 1—. (17) 

The estimates for the exponent ratios obtained are listed in Table |. The values of ^ and 
^ are within errors, for the random bond model or very close, for the anisotropic model, to 
the Ising exponent ratios f = | and ^ = |. These exponent ratios are predicted analytically 



27,28] to be of this magnitude all along the critical line of the pure Ashkin Teller model. 
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FIG. 1. The specific heat C of the random bond and the anisotropic Ashkin Teller models as 
a function of log I. The solid lines are fits to the form (|l5|). 



TABLE I. Estimates for critical exponent ratios of the Ashkin Teller models from finite size 



scaling at . obtained with linear fits to log / according to equations (15 - 17). 





a 

V 


P 

V 


1 

V 


j3W 

V 


7 (p) 

V 


fitting interval 


random bond 
anisotropic 


-0.536(32) 
-0.745(4) 


0.1252(4) 
0.1262(2) 


1.7502(5) 
1.7488(1) 


0.322(1) 
0.3312(5) 


1.3566(15) 
1.339(1) 


24 < I < 256 
16 < I < 256 



C. Calculation of Vx 

The variance a\ of the Monte Carlo estimates Xj is the sum of two contributions. The 
main contribution is due to the variance Vx of the distribution of the true X{. Vx is the 
quantity we wish to study. The second contribution is due to the errors of the estimated 



observables, 5Xi. Thus, the unbiased estimator |5T| of the variance of the X is 

Vx = a\- [(6Xl) 2 } ■ 



(18) 



5Xi depends on the length of the MC runs and on the autocorrelation time tx of the MC 
dynamics. To obtain a valid estimate of Vx, [(SXi) 2 ]/Vx should be sufficiently small. In the 
random bond Ashkin Teller model studied here this requirement was not met for the specific 
heat C, whereby we could not study Rc- Additional discussion of the practical implications 
of (Tl8| ) can be found in section III of || . Next we list results concerning the critical behavior 
of V x at T c °°. 



D. The relative variance Rx 

In Fig. [2] we plot R m , R x , Rp and R Xp as a function of log/. The solid lines are linear 
fits to the form 



R 



x 



A x l px 



(19) 



for 24 < I < 256. The estimates obtained for px are p x = —0.537(32), p m = —0.546(38), 
P x (p) = —0.493(37) and pp = —0.509(41). Clearly all observables, m, %, P, are weakly 
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FIG. 2. The relative variance of the susceptibility R x , of the magnetization R m , of the 
polarization Rp and of the polarization susceptibility R Y ( P ) of the Ashkin Teller model at as a 



function of log 10 I. The solid lines are linear fits according to (19) 
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FIG. 3. The variance of the energy Ve of the Ashkin Teller model at T£° as a function of 
log 10 L The solid line is a linear fit according to Ve = A v eI Xe , yielding xe = —2.005(26). 

self averaging. Furthermore p x and p m are in very good agreement with the value of ^ = 
—.536(32), while p Xp and pp are also within errors of 2-. Thus the results for the scaling of 
the relative variance of these four observables are in good agreement with the predictions of 
AH and with (§). 

For the energy one cannot separate the singular and the analytic parts. In addition, the 
singular part decays as l p , with p = Using our estimate of - = —.536(32) and the 

hyperscaling relation - = - — d, we find p = —1.268(48). Thus the variance of the singular 
part of the energy is expected according to @ to scale as l~ 3m2 . Therefore one would expect 
Ve to be dominated by the fluctuations of the analytic part of E decaying as l~ d . In figure 
[3] the variance of the energy Ve as a function of log / is plotted. Straight forward linear fits 
to the form Ve = A v eI Xe in the lattice size range 24 < / < 256 yielded xe = —2.005(26) in 
good agreement with our expectation xe = —d. 

To conclude this part of the study, we found weak self averaging at criticality for a 
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disordered model governed by a pure fixed point with (^j < 0. We also found good 
agreement with the scaling prediction (0). 

III. THE SITE DILUTE ISING MODEL ON A CUBIC LATTICE 

The second model we chose to study is the site-dilute Ising model on a cubic lattice ( 
see e.g. and references therein). On every site of a I x I x I cubic lattice either an Ising 
magnetic spin Si — ±1 is placed if Ki — 1 or a vacancy is placed if Ki — 0. The Ki are 
randomly drawn according to one of the prescriptions given below. The system is governed 
by the Hamiltonian 

U = - J ]T KiSiKjSj , (20) 

<i,j> 

where < i,j > stands for a pair of nearest neighbors. RG calculations found a dilution 
independent random fixed point with universal critical exponents. For example, a recent 
calculation @ obtained 7 = 1.313, (3 = 0.342 v = 0.666 and a = 0.002. Early MC 



studies found global effective critical exponents which were found to depend on dilution. 
This was later interpreted as due to crossover effects. For example, in a most extensive 
MC study, Heuer []32| found from finite size scaling in the lattice size range 20 < / < 60 
values ranging from -{p = 0.95) = 0.12(6) with 5% vacancies, to -(p = 0.8) = —0.04(6) 
and ~(jp = 0.6) = —0.22(6). However he argued by analyzing a suitable scaling function, 
that all models with different amounts of dilution are exhibiting a crossover to the fixed 
point predicted by RG. His results show that, of the amounts of dilution he studied, the 
p = 0.8 model reached the universal behavior at the smallest lattice sizes. Later Janssen 
Oerding and Sengespeick ]20| showed in their RG calculations that the effective exponent 
values obtained by Heuer can be related to regions in the space of coupling coefficients away 
from the fixed point. 



A. Details of the simulations 

Three site-dilute Ising models were examined, including two types of disorder. In one 
model disorder was realized in a canonical manner; namely, the number of magnetic sites in 
each site-diluted sample was fixed at a fraction c = 0.6 of the number of sites in the lattice. 
Thus fluctuations among samples occur only in the locations of the magnetic sites but not 
in their number. In two other models disorder was realized in a grand-canonical manner; 
namely, each sample was created by assigning to each site of the lattice a magnetic spin 
(vacancy) with probability p (1 — p). In one model we used p = 0.6 and in the second one 
p = 0.8. In this case fluctuations among samples include fluctuations in the number of mag- 
netic sites. These fluctuations tend to zero as I — > 00, but for finite I they are significant. For 
this reason we found it of interest, in this study of fluctuations among samples, to compare 
the two ensembles. We are unaware of any previous findings attesting to differences in the 
asymptotic critical behavior between the two ensembles. Because of the (spatially) uncor- 
related nature of the disorder in the grand canonical ensemble it is favored for its relative 



simplicity by theoretical studies (see for references) and by numerical studies [|33| -|37 
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TABLE II. Number of site-dilute samples simulated for each model and lattice size I. The 



last column lists the infinite lattice critical temperature, as estimated by Heuer [32|, at which the 



simulations were performed. For c = 0.6, I = 90 the pseudocritical temperature was not estimated. 







I = 10 


I = 20 


I = 40 


I = 60 


I = 80 


I = 90 


1 C 


c = 


0.6 


8000 


26000 


2000 


800 




1000 


2.4220(6) 


p = 


0.6 


72000 


47000 


8000 


950 


800 




2.4220(6) 






I = 4 


I = 8 


I = 16 


/ = 32 


/ = 64 






p = 


0.8 


10000 


4000 


32000 


4000 


1479 




3.4992(5) 



aiming to test them. On the other hand, in studying by Monte Carlo average thermody- 
namic observables, errors can be reduced by using canonical disorder, as was done in [R2j . 



We note that if one wishes to study by Monte Carlo the fluctuations in the thermodynamic 
observables due to disorder, the use of grand-canonical disorder is advantageous. 



In the Monte Carlo simulations we used the Wolff [B(J single cluster algorithm [29] for the 



Ising model because of its efficiency [B8J . Skewed periodic boundary conditions M were used 
in order to speed up the simulations. For each model and lattice size n site-dilute samples 
were simulated. Table [I| summarizes the number of samples n used for each lattice size 
for the three models. Simulations were performed at the estimated infinite lattice critical 
temperatures T c °°, given in table |J due to Heuer [fj^] , and taking J = 1. The procedure for 
calculating the T c (i,l) is described in Sec. |V A . 



B. Finite size scaling at T c °° 

Here we give the finite size dependence of averages (over all samples) of some observables 
at the critical point [Xi(T£°)]. 

1. Magnetization m and susceptibility \c 

Using the abbreviation 

M = Y,K i S i , (21) 

i 

the magnetization density m is defined as 

m =», ,22, 

where N = I 3 and the fraction of magnetic sites was either p = 0.8 or p = 0.6. The 
susceptibility at T c °° was denned as 

V - ^ (231 
The magnetization density m and susceptibility Xc were fitted to the finite size scaling forms 
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form fll6|), yielding estimates for — which are listed in table III 
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TABLE III. Estimates for order parameter critical exponent ratios from finite size scaling at 
T^°. Estimates due to Heuer are listed for comparison. Note that throughout the paper the 
errors given for our results are only statistical while Heuer's error estimates include the systematic 
errors which could arise from errors in determining T^°. 





§_ 

V 


f; (Heuer) 


1 

V 


^ (Heuer) 


c = 0.6 


0.438(13) 


0.45(2) 


2.110(21) 


2.09(3) 


p = 0.6 


0.437(12) 




2.104(20) 




p = 0.8 


0.505(2) 


0.51(2) 


1.990(4) 


1.98(3) 



The estimates which were obtained for the critical exponents ratios - and - from the 
fits in figures f| and [| are listed in Table |TTT| together with the estimates of Heuer [32|| . Note 
that exponent ratios and critical temperatures quoted from [32] for p = 0.8 were actually 
obtained with canonical disorder c = 0.8. 



2. Tjjr and estimation of ^ 

In attempting to estimate the exponent ratio - directly from the finite size scaling of 
the specific heat C through (|l~5|), we encountered two difficulties. First, we found that the 
estimates of the specific heat of each sample were very sensitive to the length of the 
simulations. Shorter simulations biased the specific heat to lower values, e.g. for p = 0.6, 
I = 40 The value measured for C using average simulation length of nj (2te + 1) ~ 120 was 
two standard deviations smaller than the one measured with a four times longer simulation. 
The systematic underestimation of response functions due to run lengths which are too short 



was studied in f39|. Second, the accuracy in estimating ^ from the specific heat behavior is 
rather poor. This is due to the fact that ^ is a small negative number so that the singular 
behavior of C is difficult to disentangle from other analytic contributions ^ . 



In order to overcome the difficulty in estimating the exponent ratio - we followed Heuer 
40| and measured the derivative of the magnetization with respect to the reduced temper- 



ature t = ?-7f^. It is equal to the magnetization-energy correlation 

J- c 

r = -7ST = -Wm *> Wr* [{m ~ < |M|>)(H " {n)m • (24) 

where / is the free energy density. From the scaling behavior of the free energy 

f{t,h) = b- d f{W%by*h) (25) 
one finds that Y diverges as t~^, where 

( = yt - {d - yh) =1-P. (26) 

Thus we fit T to the finite size scaling form 

T = C,li . (27) 
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TABLE IV. Estimates for critical exponent ratios from finite size scaling at T£°. Estimates due 
to Heuer |3^] are listed for comparison. The estimate for ^ is based on the relation ^ = 2(£ + ^) — d. 







£ (Heuer) 


(?) 


2 (Heuer) 


c = 0.6 


0.948(6) 


0.94(2) 


-0.228(28) 


-0.22(6) 


p = 0.6 


0.958(3) 




-0.210(30) 




p = 0.8 


0.962(4) 


0.97(2) 


-0.066(9) 


-0.04(6) 



The resulting estimates for £ are given in table 0. Assuming the hyperscaling relation 
2 = 2 — d and using ( p6|) the scaling relation ^ = 2(£ + ^) — c? is obtained. Using this 
relation, the results for m and V are utilized in table [IV] to give estimates for ^ which are 
much more accurate than those obtained from analysis of the specific heat results. 



IV. LACK OF SELF AVERAGING AT T«2F 



In order to obtain the variance Vx and the relative variance Rx the same procedure and 



considerations as described in section [II C| were used. 

In figure |6| we plot the relative variance of the magnetization R m as a function of lattice 
size on a double-logarithmic scale. Several interesting features are suggested by this figure. 
First, note that for p = 0.6, R m is decreasing as I increases for the smaller lattice sizes, 
possibly leveling off for large I. R m of the p = 0.8 model first decreases slightly and then 
seems to tend to a constant. Since it seems plausible that R m {p = 0.6) > R m {p = 0.8) 
for any lattice size, these trends seem to imply that for the two grand-canonical models 
R m tends to the same constant. Assuming that this constant is bound from above by the 
p = 0.6 model and from below by the p = 0.8 model we estimate it as R m = 0.055(2). 
The implication of this scenario is that R m of the weakly diluted p = 0.8 model reaches the 
universal R m value of the dilute Ising fixed point at smaller system sizes than the highly 
diluted p = 0.6 model. The fluctuations in m; in the highly diluted p = 0.6 model are larger 
than those of the dilute Ising fixed point model. This finding is in line with Monte Carlo 



results |32| and RG calculations |20| , according to which the critical exponents of the dilute 
Ising fixed point are closer, in the lattice size range 20 < / < 60, to the observed effective 
critical exponents of the p = 0.8 model than to those of the p = 0.6 model. 

A second feature is the striking difference between the two types of disorder, with the 
canonically disordered c = 0.6 model exhibiting a much smaller relative variance than that 
of the two grand-canonical disorder models. While R m of the c = 0.6 model is initially 
increasing with system size it appears to level off to a constant value of R m {l = 90) = 
0.0227(8). Though it is possible that R m could increase at larger lattice sizes it seems 
unlikely since the system sizes are already quite large. An indicator to the similarity of the 
two types of models p = 0.6 and c = 0.6 is the relative square root mean of the fluctuations 

in the number of magnetic sites M = J2jLi Kj in the p = 0.6 systems \J ^y$£^ — \j 

which for I = 80 is as small as ~ y^j. 

If indeed R m of the c = 0.6 model tends to a different constant than that of the models 
with grand canonical disorder, then according to Aharony and Harris' very general RG 



arguments Jl3| the two types of models do not belong to the same universality class! We 
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FIG. 6. The relative variance of the magnetization R m at as a function of log 10 ^- 



are not aware of any additional evidence to this effect otherwise. For example our critical 
exponent estimates for the c = 0.6 and p = 0.6 are compatible with each other, and our 
exponents for the p = 0.8 model are compatible with those of Heuer [32] for a c = 0.8 model. 
The critical temperatures for both types of models seem also to agree (see table [V] and 
References [32,35|)- This question is currently under investigation. Preliminary results fr4lfl 
suggest that the two types of models do flow to the same fixed point and that the difference 
in R m will disappear for very large I. 

The relative variance of the susceptibility R Xc is plotted in Fig. [T[ R Xc exhibits the same 
qualitative behavior as that of R m . R Xc of the grand canonical disorder models tends to 
R Xc = 0.156(4), while R Xc of the c = 0.6 model seems to tend to R Xc (l = 90) = 0.061(2). 
Aharony and Harris |13|J42]| found that to leading order in e = 4 — d, Rm/R x = 1/4. We 
find that for p = 0.8 R M /R X = 0.35(2) and for c = 0.6 R M /R X = 0.37(1). Possibly terms of 
higher order in e would reconcile this discrepancy. It cannot be attributed to the definition 
of Xc- If one defines the susceptibility as in ( ffij| ) then at T c °° one finds that R x becomes 
smaller by a factor of ~ 7 — 10. In this case the ratio Rm/R x would become even larger. 
We did not use this definition for the susceptibility at T c °° because of its large single sample 
errors 8xi ( see a l so PI)- 



V. SCALING OF PSEUDO CRITICAL TEMPERATURES 

A. Calculating T c (i,l) with the histogram reweighting method 

One of the main purposes of this work was to study the distribution of pseudo-critical 
temperatures T c (i, I) of the ensemble of site-dilute Ising models. The main aim was to 
study directly the scaling of 5T C {1) with / and test which one of (j5[) or (|8|) is correct in 
the case of a system governed by a disordered fixed point. The inverse pseudo-critical 
temperature K c (i,l) = 1/T c (i,l) of the z'th sample was defined as the inverse temperature 
of the maximum of the susceptibility of that sample, K c (i, I) = K max (i, I). Here the definition 
of the susceptibility was 
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FIG. 7. The relative variance of the susceptibility R Xc at T c °° as a function of log 10 I. 

(M 2 ) - (|M|) 2 



Xi 



NpT 



(28) 



In order to find K max (i, I) the following iteration procedure was followed for each sample. 
A first simulation was performed at the infinite lattice critical temperature (as estimated 
in [32]) Ki = Kf (the index i is omitted from here on). In addition to calculating the 
observables m, x, T, a histogram of the energy and magnetization was generated. Using the 
single histogram reweighting technique pl|-p5H (For previous studies of disordered systems 
utilizing the histogram reweighting technique see PB] , ^5[ ), this histogram can serve to calcu- 
late observables at temperatures close to K\. By calculating \ at different temperatures a 
first estimate for the susceptibility maximum x m ax an d the temperature at which it occurs 
-^max was obtained. A second simulation was then performed at a temperature somewhat 
above this estimate K 2 = K^ ax — of f set. Previous studies of the histogram reweighting 
technique have shown that the errors of observables at T, 5X(T), are smaller [^,^3|] when 
the temperature at which the histogram was generated T sim is slightly higher, T sim > T. For 
this reason we chose offset to be a small positive number (see more details below). 

Using the energy and magnetization histogram generated at K 2 a new estimate for the 
temperature of the susceptibility maximum, ^ ax , was obtained. If the difference between 
the two estimates was smaller than a predetermined resolution r, \K^ ax — K^ ax \ < r, the 
iteration process was stopped. Otherwise the iteration process continued, whereby Kj = 
-^max ~ °f f se ty until the condition 



\Ki -TP™ 1 ! 

I max max I 



< r 



(29) 



was met. This iteration process was intended to overcome the problem of systematic errors 
23| that occur when the simulation temperature is too far from the true K maiX . The condition 
( P^D is supposed to ensure that the last two estimates for K max do not suffer from a systematic 
error, r was chosen equal to the approximately expected statistical error of i^4ax- ^ the 
iteration process did not terminate before or with the third estimate K^ ax then the Monte 
Carlo simulation length at was doubled and the process was continued. It was again 
doubled if it reached the seventh iteration and again doubled if it reached the tenth iteration. 
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Non-convergence of the process after twelve iterations was very rare. In those samples the 
iteration procedure was restarted manually with K\ = but with a larger initial Monte 
Carlo simulation length. The need to increase the simulation length for some samples occured 
because for different samples there were different auto-correlation times (of the Monte Carlo 
dynamics) and different average cluster sizes, while the simulation length was specified by 
the number of Wolff cluster flips. 

In order to estimate the statistical error and reduce it, a simulation with five times 
as many Monte Carlo steps (compared to the simulation length of the last iteration) was 
performed again at the last simulation temperature Kj. The Monte Carlo sequence was 
broken into five, using each segment to create a separate histogram and calculate a separate 
estimate of K mauX and Xmax- Together with the last estimates of K max and Xmax of the 
iteration procedure, all together six estimates of K m3jX and Xmax were averaged to give final 
estimates of K maXi i and Xma X ,«- The variance of these six estimates was used to estimate the 
error for the two quantities, 6K maXi i, and ^Xmax,*- 

The parameter offset was adjusted for the small system sizes, through trial runs, so 
as to minimize the errors in Xmax, while its value for the larger systems was extrapolated 
from the smaller ones. For c = p = 0.6 we set offset « 0.27Z~ L66 , and for p = 0.8 
offset « 0.12/ -1,63 . The optimal value of offset was found not to depend strongly on the 
simulation length. The resolution r was adjusted so as to be approximately equal to the 
ensemble average statistical error of -K^ ax . Note that the parameters offset and r were set 
once for each model and each lattice size and were not varied for different samples. In some 
of the larger systems, for a subset of the samples, the simulation with five times as many 
Monte Carlo steps was not performed, so that error estimates of K max and Xmax were not 
obtained. This was done in order to save computer time. For these samples the average 
squared error of K max and x max was approximated as being six times larger than that of the 
complementary subset of samples where the error was calculated (from an all together six 
times longer Monte Carlo sequence). For the p = 0.6 I = 80 system the estimated average 
squared error was extrapolated from the smaller systems. 



B. Scaling of t c (l) 

In the finite size scaling theory of |9j it was assumed that the average pseudo-critical 
temperature K max (l) = [K max (i, I)] scaled as 

^max(0 ~K C = A K r X , (30) 

and that the shift exponent |4|| A = y t = 1/v. First we assumed the correctness of the 
critical temperature values T c °° quoted in table [0], so that the critical inverse temperature 
of the infinite sample is assumed to be K c = 1/T C °° = .285781(40) for p = 0.8, and K c = 
0.41288(10) for c = 0.6 and p = 0.6. Fitting K max (l) to © with K c fixed we found for the 
p = 0.8 model values of A = 1.446(34) and A K = 0.040(5) from lattice sizes 16 < I < 64. 
For the c = 0.6 and p = 0.6 models the results were incompatible with the fixed value of 
K c = 0.41288. In fact in these models K max (l) monotonically decreases with / and for the 
largest lattices we have K max (80,p = 0.6) — K c = —0.00026(3) and i^ max (60, c = 0.6) — K c = 
—0.000077(35). Thus we also fitted i^ m ax(0 to (PH|) with K c being a free parameter. The 
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TABLE V. Different parameters related to the average pseudo-critical inverse temperature 
K max (l) and its variance Vjr max . Second column: estimate of the shift exponent A according to 



pOl ), where K c is taken from Heuer [32]. Third and fourth column: same as first column but with 



K c being a free parameter. Fifth column: estimate of yt based on the finite size scaling of m and 



r. Sixth column: same as fifth according to [32]. Last column: exponent of Vk, 





A, (K c fixed) 


A 


K c 


yt = i + ? 


y t (Heuer) 


PK 
2 


c = 0.6 




0.99(19) 


0.41254(13) 


1.386(14) 


1.39(4) 


1.41(4) 


p = 0.6 




1.30(10) 


0.41251(5) 


1.395(12) 




1.421(9) 


p = 0.8 


1.446(34) 


1.346(2) 


0.2857609(4) 


1.467(5) 


1.47(4) 


1.44(2) 



values of A and K c which were found, using lattice sizes 10 < I < 60 for c = 0.6 and p = 0.6, 
and 8 < / < 64 for p = 0.8, are given in the third and fourth columns of table |V|. For 
p = 0.8 our estimate K c = 0.2857609(4) is within errors of the estimate of Heuer |32] (with 
canonical disorder) and of Wang et al |3S] (with grand canonical disorder). For c = 0.6 and 
p = 0.6 our estimates K c = 0.41254(13) and K c = 0.41251(5) are within errors of each other 
but not within errors of the assumed value K c = 0.41288(10). A more accurate estimate of 
K c , which does not require knowledge of A, is obtained in section |V C 1| . In the fifth column 
of table [V] estimates of y t based on estimates of -, - and the scaling relation (|26|) are given. 
The values of y t obtained by Heuer in the same way are given in the sixth column of table 
[V] and are compatible with our estimates. Our estimates of y t and A agree for p = 0.8 
(where K c was fixed) and for p = 0.6 (where K c was a free parameter). For c = 0.6 no 
agreement was found. One possible reason could be that the system size used to estimate 
A was too small and that corrections to scaling need to be taken into account. As is well 
known, finding the critical temperature and the shift exponent simultaneously is a difficult 
task. In any case, our estimates for y t are much more accurate than the estimates for A. It 
seems to us that trying to extract the shift exponent and the critical temperature by finding 
the pseudo-critical temperature of many samples and using their average K max (l) in (|30|) is 
not an efficient method. This is because a long MC simulation is needed to avoid systematic 
errors in estimating K max (i,l). Thus it is difficult to obtain a sufficient number of samples 
for an accurate enough estimate of K max {l). The error in K max (l) must be small compared 
to K max (l) — K c . This difficulty is not as significant for the estimate of the variance Vft; max . 



C. Variance of pseudo-critical temperatures V/f max 

The variance of the pseudo-critical temperatures distribution Vx max was calculated taking 
into account the errors [(SK maXti ) 2 }. This is done in a manner completely analogous to the 
discussion of Vx in Sec. |II Q . Vx max is plotted in Fig. || on a double logarithmic scale. The 
solid lines are fits to the form V^ max = A^l~ pK and the resulting estimates of ^ are listed 
in the last column of table |V]. As one would expect, V# max is smaller for c = 0.6 than for 
p = 0.6, and is the smallest for p = 0.8. We see that for all three models the results for px 
exclude the possibility (|5|) that px = d = 3. On the other hand 4p is within errors of yt for 
p = c = 0.6, and within errors of A (with K c fixed) for p = 0.8, as predicted by Aharony and 
Harris (§). Note that the values obtained for p K for p = 0.8 with lattices 8 < / < 32 and 
p = 0.6 with 20 < / < 60 are p K = 2.95(6) and p K = 3.00(4). This behavior of V Kmm could 
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FIG. 8. The variance of the inverse pseudo-critical temperatures V#- max as a function of I on a 
log — log scale. The lines are linear fits yielding exponents px listed in the last column of table |V|. 

be a manifestation of a crossover from pure @ to dilute (H) critical behavior. On the other 
hand for the model with the canonical disorder the crossover is in the opposite direction 
since for c = 0.6 with 10 < I < 40, p K = 2.77(7). 

The results for V/<: max support the picture implied by AH RG calculations, namely that 
both the width of the pseudo-critical temperatures JVk^JT) and the distance of its average 
from the critical inverse coupling \K max (l) — K c \ scale as ~ l~ yt . This is best visualized in 
Fig. H where the frequency of the scaled pseudo-critical inverse temperatures (K max (i,l) — 
K c )l yt is plotted for p = 0.8 and I = 16,32,64 with K c = 0.285781 and y t = 1.467. It is 
evident that the three distributions match well. Their averages are [(K max (i,l) — Kc)l yt ] = 
0.047(1), 0.0451(28), 0.044(5), and their widths are \/7^l yt = 0.172(15), 0.174(28), 0.17(4) 
for / = 16, 32, 64 respectively. Note that the average ratio of the width to the average is 
ps 3.8 . Thus, as is evident from Fig. [5|, the fluctuations in K c (i, I) are significantly larger 
than \K c (l) — K c \ for any system size. The result is that a measurement of X at the critical 
temperature T c °° is done in some samples above their pseudo-critical temperature T c (i, I) 
and in some samples below T c (i, I) ! 



1. Estimating T c °° through Vk mtre 

Our estimates of V^ max allow us to estimate T c °° by another method (We thank D. Stauffer 
for bringing this to our attention). Since asymptotically K milx (l) — K c ~ l~ yt and y/Vjc^ ~ 
l~ yt , one expects that 

K mSK (l) = K e + B v y/V Knm (l), (31) 



where K c and B v need to be determined. Note that by fitting the data according to (JJT) 
(this method was used in percolation studies |45]] ) it is not necessary to determine z/, and 
only two fitting parameters are used. Therefore the estimates of K c obtained in this way are 
probably more reliable than those given in table |V[ In figure [10] we plot K mSuX as a function 



°f V^Kma* together with linear fits made according to (pTf) . We find K c = 0.285779(2), 
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B v = 0.257(3) for p = 0.8, K c = 0.41244(4), B v = 0.22(3) for p = 0.6, and K c = 0.41265(4), 
B v = 0.230(15) for c = 0.6. 




D. Maximum of the susceptibility Xmax 
1. Scaling of [x 

Another way to study the finite size scaling of the susceptibility is to study the ensemble 
average of the maximum susceptibility [Xmax] which is expected to scale with lattice size 
as in (0) with a scaling exponent ^. In figure [Xmax] is plotted as a function of I on a 
double logarithmic scale. The straight lines are linear fits to the form (jIB[). For p = 0.8 we 
find - = 1.987(3) which is in agreement with the estimate obtained using the susceptibility 
Xc at T 6 °°, 1 = 1.990(4). As can be seen in Fig. |TT] the values of [x max ] for V = 0.6 and 
c = 0.6 are indistinguishable (they are indeed within errors). This is in contrast with the 
data at T c °° of fig. |5| where Xc of the two models seem to diverge with a similar exponent 
but with a different amplitude. We have also calculated x a t and found the same trend, 
namely that x{p = 0-6) > xi c = 0-6), so that this feature is not an artifact of the different 
definitions, (|23|) and (p8j), for x- F° r V — 0.6 and c = 0.6 we found ^ = 2.027(2) and 
- = 2.034(2) respectively. These values are significantly lower than the values found using 
the susceptibility Xc at T c °°, ^ = 2.104(20) and ^ = 2.110(21). They are also closer to results 



of RG calculations £ = 1.97 fl20|JT 
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FIG. 11. The ensemble average of the maximum susceptibility [% n 
The solid lines are fits to the form (jltf), yielding estimates of ^ = 2.034(2) for c 
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2. Lack of self averaging of the relative variance 



In figure [T2| We plot the relative variance of the maximal susceptibility -R Xmax as a function 



of lattice size on a double-logarithmic scale. For p = 0.8, I = 64 and p = 0.6, I = 80 we 
have [(<5xmax,i) 2 ]/^:max ~ | , 0.15 respectively (see section [II Q ). Thus the estimate of V Xm ^ 



is dominated by the estimate of the average squared single sample errors [(<5x max ,i) 2 ]- F° r 
p = 0.6 I = 80 this estimate was actually extrapolated from the smaller systems estimates 
(see | V A| ) . Thus these two data points should be taken with more than a grain of salt. It is 
most interesting to compare Fig. [12| to Fig. [7] where the relative variance of the susceptibility 



at T c °°, R Xc , is plotted. For p = 0.8 the behavior of -R Xmax and R Xc is qualitatively rather 
similar with -R Xmax initially decreasing as I increases and tending for larger I to a constant, 
where -R Xmax (^ = 64) = .00216(16). However, in contrast with Fig. [7], this constant is roughly 
72 times smaller than the large I value of R Xc . This is quite a striking difference. It means 
that in order to obtain the same relative accuracy in [% c ] as in [x max ] approximately 70 
times as many samples are needed. The source of this difference is apparently simple. The 
susceptibility of each sample is some function G of the temperature with a sharply peaked 
maximum at T c (i, I). In fact G is approximately only a function of the difference T — T c (i, I), 
G(T — T c (i, I)). Thus, the value of the maximum susceptibility is nearly sample independent, 
Xmax ~ G(0). On the other hand, when one measures x a ^ T£°, in different samples one is 
sampling G at different values of its argument. This results in large fluctuations in x at T c °°. 

Our findings suggest that the standard procedure of investing much computation time 
in finding the I —>■ oo limit of the critical temperature, T c °°, and then averaging quantities 
at this temperature over many samples, is not optimal. A better procedure may be to 
locate through the single or multiple histogram method the pseudocritical temperature of 
each sample, and measure quantities at that temperature. In this way sample to sample 
fluctuations are reduced substantially and better accuracy is achieved. 

For p = 0.6 -R Xmax monotonically decreases with lattice size, possibly leveling off to a 
constant for large In contrast with R Xc , this constant is different from that of the p = 0.8 
model. Lastly, for c = 0.6, in contrast with R Xc , we find that -R Xmax is within errors of -R Xmax 
of the p = 0.6 model. In addition, R Xmax initially decreases as / increases, opposite to the 
behavior of R Xc . More explanations to the differences in the behavior of -R Xmax and R Xc are 
given at the end of the next subsection. 



E. Dependence of m(T c °°) on T c {i,l) 

After examining the behavior of the distribution of X(i, I) at T c °° and the distribution 
of T c (i, I) it is imperative to examine the correlation between the two distributions. A 
good starting point is the finite size scaling ansatz (|j), according to which Xj(T c °°) mainly 
depends on Tc Fig. [13| is a scatter plot where for each sample % the horizontal axis 



represents the scaled absolute inverse temperature \K C — K c (i,l)\l yt and the vertical axis is 
the scaled magnetization m^" . This representation is equivalent to the usual data collapse 
representation which is used to demonstrate finite size scaling. The difference is that here the 
reference critical temperature is K c (i, I) instead of K c , and the measurement temperature is 
always K c instead of different values of K. Points with K c > K c (i, I) constitute the higher 
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FIG. 12. The relative variance of the maximum susceptibility -R Xmax as a function of log 10 I. 



m (lower temperature) branch, whereas points with K c < K c (i, I) constitute the lower m 
(higher temperature) branch. In figure |13| we plot data for p = 0.8 and I = 16, 64. For 
the sake of clarity, only 100 points are shown for each system size and each branch, and 



several points with \K C — K c (i,l)\l yt < 0.001 were omitted. Figure [13] indicates that to a 
good approximation the scaled magnetization of the sample at T c °° is a function of only the 
scaled reduced temperature of the sample. Thus one may attempt to substitute @ by a 
sample independent form for Qi(Z) so that 

Xi(T,0~Z p g(*iP) ■ (32) 

Note that this is only a good approximation; if ( j32|) were exact, it would mean that -R Xmax = 
0. Thus in order to describe the magnetization data at K c we write (the change from 
temperature to inverse temperature is only for convenience) 



mi (K,l) =r-Q ± {\(K - K c (i,l))\l*} 



(33) 



Here Q + (Z) is the scaling function for K < K c (i,l) and Q_(Z) for K > K c (i,l). For 
large I, and thus large Z, the infinite sample critical behavior, ~ {K — K c (i,l)} , must 
be asymptotically reproduced for K > K c (i,l). Thus, for large Z, Q-(Z) ~ Z 13 . 

decay of the 

K c (i,l), 



For K < K c (i,l) the shape of Q + (Z) must reproduce, for large Z, the 



/N 

~ ft - 

magnetization to zero as I — ► oo. Thus, for large Z, Q + {Z) ~ Z 2yt . For K 
i.e. Z = 0, the finite size scaling behavior [m*] ~ l~» must be asymptotically reproduced^, 
implying Q±(Z) — > const as Z —>■ 0. A simple possible form for Q±(Z) fulfilling these 
requirements is 



Q±(Z) = A±Z P± {1 + B±z~ p± y± 



(34) 



where p_ = (3 and p + = (3 - 
13| should be described by 



2^ and A±, B±,p± are free parameters, so that the data of Fig. 
J34]) with = fT c . Thus, for each lattice size I = 16, 32, 64 



An example of this type of finite size scaling is the scaling we found, for 
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FIG. 13. Scatter plot, where for each sample i the horizontal axis represents the scaled 
absolute inverse temperature Z = \K C — K c (i, l)\l Vt and the vertical axis is the scaled magnetization 
Q± = ml" . Points with K c > K c (i, I) constitute the higher m (lower temperature) branch, whereas 
points with K c < K c (i,l) constitute the lower m (higher temperature) branch. For the sake of 
clarity, only 100 points are shown for each system size and each branch, and several points with 
\K C - K c (i, l)\P* < 0.001 were omitted. 

TABLE VI. parameters of the fitting functions Q±(Z), defined in (pi)), obtained by fitting the 
data sets of scaled {K max (i, I), rrii(K c )} pairs for each lattice size I = 16,32,64 separately (with 
p=0.8). The number of samples used was 32000, 4000, and 1477 for I = 16, 32, 64 respectively. The 
six fitting functions are plotted in Fig. [l4|. The crossover lengths, which control the crossover to 



the large Z behavior are defined as Z± oss = B± 





A_ 




P- 


^•cross 


A + 


B+ 


P+ 


^cross 


I = 16 
I = 32 
I = 64 


2.387(4) 
2.380(9) 
2.291(23) 


0.0518(13) 

0.048(3) 

0.07(1) 


1.45(1) 
1.49(3) 
1.33(6) 


0.130 
0.130 
0.142 


0.4695(12) 
0.4623(22) 
0.4524(75) 


0.1765(18) 
0.1687(34) 
0.152(11) 


1.299(5) 
1.316(9) 
1.368(35) 


0.263 
0.259 
0.252 



(a partial set of which is plotted in figure [13]) were fitted to the form (p4|). The values of 
p_ = (3 = 0.34295 and p + = (3 — ^ = —0.67565 which were used rely on the finite size 
scaling results at T c °° (tables |T| and 0). The six fitting functions which were obtained 
are plotted in Fig. [H] and their fitting parameters are given in Table 0. The agreement 
between the three curves for both branches, as seen in Fig. O, is surprisingly good. The 



and both branches, K c < K c (i,l) and K c > K c (i,l), the scaled {K max (i, I), mi(K c )} pairs 



goodness of the fits is also extremely high. This suggests that equation ( p2[ ) , equations 
similar to fl3~4l) , and the possibly invariant (as suggested by Fig. || ) distributions of K c (i, I) 
provide an excellent description of the scaling behavior of disordered systems. 

In figure 15 we show the same data as in figure 13 but for p = 0.6 and c = 0.6 with system 
size I = 40. The purpose of this analysis is to demonstrate that the magnetization of the two 
models is governed by the same temperature dependence, and that the main difference is in 
the distributions of K c (i, I). For this reason the data were scaled with the same exponents, 
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FIG. 14. The functions Q±(Z), as defined in (]3"4|), obtained from best fits to the scaled 
magnetization versus temperature scatter plots for I = 16, 32, 64. upper curves according to Q- 
(K c > K c (i, I)) and lower curves according to (K c < K c (i, I)). The fitting parameters are given 
in table [Vj 
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FIG. 15. Same data as in figure |13| but for p = 0.6 and c = 0.6 with system size I = 40. 
The data were scaled with exponents taken as the average of the exponents of the two models, 
^ = 0.4375 and yt = 1.3905. For the sake of clarity, only a 100 points for each model and each 
branch are shown. 
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0.001 0.010 0.100 1.000 

FIG. 16. The functions Q±(Z), as defined in (|34|), obtained from best fits to the scaled 
magnetization verses temperature scatter plots for c = 0.6 (dotted line) and p = 0.6, with I = 40. 
upper curves according to Q- (K c > K c (i,l)) and lower curves according to Q + (K c < K c (i,l)). 
Fits made using ^ = 0.4375 and yt = 1.3905. 



taken as the average of the exponents of the two models, ^ = 0.4375 and y t = 1.3905. In 
fact our estimates for y t and ^ for the two models are within errors. For the sake of clarity, 
only 100 points for each model and each branch are shown. As was seen with the p = 0.8 
data, it is evident that to a good approximation in both models the magnetization at K c is a 
function of only the reduced inverse temperature K c — K c (i, I). The main difference between 
the two models is also clear; For p = 0.6 there are more points with large \K C — K c (i,l)\, 
while for c = 0.6 there are more points with small \K C — K c (i, Thus larger fluctuations 
for p = 0.6 in K c (i,l) (see also Fig. §) together with the large dependence of rrii(K c ) on 
K c — K c (i, I) give rise to the result that R m (p = 0.6) > R m (c = 0.6). 

In figure [16] we plot the fitting functions Q±(Z), obtained by best fits to the scaled 
magnetization verses temperature scatter plots for p = 0.6 and c = 0.6 with I = 40 (the full 



data sets corresponding to figure [T^). For the high temperature branch (K c < K c (i, I), lower 
curve) good agreement between the fitting functions Q+(Z) of the two models is found. For 
the low temperature branch (K c > K c (i, I), higher curve) good agreement is found between 
the functions Q±(Z) for smaller Z, while for large Z Q_(Z) is larger for the grand canonical 
disorder (p=0.6). The fitting functions Q±(Z) for the data for I = 60 did not agree with 
those of I = 40. Possibly this is so because the exponents used are not the asymptotic ones 



It is also of interest to contrast the dependence of Xmax on K maxi with the dependence 
of Xc{K c ) on K max j. This is a key to understanding the reasons for the differences between 
the characteristics of R Xc (figure 0) and the characteristics of -R Xmax (figure [12]) . In figure |I7| 



we show a scatter plot of (-Kmax — K c , ^ max j ) and (K max — K c , t^|Ht) f° r V = 0-6 and system 
size I = 60 from 950 samples. It is evident that while Xc(K c ) shows a strong dependence on 
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FIG. 17. A scatter plot of (K max - K c , ) and (K max - K e ,gj^), contrasting the 

dependence of Xc(K c ) and Xmax on K max (i,l). Data for p = 0.6 and system size I = 60 from 950 
samples. 

K max — K c , Xmax shows little dependence on K max — K c . This qualitative difference persists 
for all models and all system sizes. This explains why, for any given model, fluctuations 
in i^max.i give rise to fluctuations in Xc{K c ) which are much larger than the fluctuations in 
Xmax- The result is that i? Xmax -C R Xc , as we have noted previously. 

Fig. [17] is also the key to understanding why R Xc (p = 0.6) > R Xc (c = 0.6) while 
i? Xmax (p = 0.6) ~ -R X max( c = 0-6)- In the first case, since fluctuations in i^ max are larger for 
p = 0.6 than for c = 0.6 (see figure |8|) the strong dependence of Xc(K c ) on K max — K c gives 
rise to R Xc {p = 0.6) > R Xc (c = 0.6). In the second case, despite the fact that fluctuations 
in .Kmax are larger for p = 0.6 than for c = 0.6, the weak dependence of Xmax on -ft" max — K c 
results in i? Xmax (p = 0.6) w # Xmax (c = 0.6). 

VI. SUMMARY AND DISCUSSION 

By and large it seems that our MC results confirm the AH scenario. In an Ashkin- Teller 
model, governed by a pure fixed point, we found that Rx ~ l~ ■ in agreement with (§) and 
(0). In site dilute Ising models on a cubic lattice, governed by a random fixed point, we 
found a lack of self averaging for both canonical and grand canonical disorder. One of the 
aims of our work was to resolve whether at random fixed points our assumption (|5|), which 
led to the prediction (|) for the critical width Rx, is correct? The alternative Rx — > const 
result of AH implies that (§) should replace (|). Our results indicate that the AH result is 
the correct one. Note though that the absolute value of the exponent ratio - of the dilute 
Ising fixed point, either as calculated by RG, ^ = 0.003, or as indicated by the p = 0.8 
results ^ = —0.055(8), is very small. Thus one could argue that our results for R m and R Xc 
do not disprove (|]). The scaling of VR max is, however, in agreement with @ and not with 
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(§). This therefore rules out (|J) since it is based on (|). 

We find it appropriate to repeat here the results of Aharony and Harris [13], which we 
have now validated, with an emphasis on the implication to experiments. In finite size 
scaling form the relative variance can be written as 

R x (U) = l w Q(l/i) ■ (35) 

For a fixed £ = £ and / >> £, and thus large Z, strong self averaging, R x (£ ,l) ~ l~ d , 
must be asymptotically reproduced. Thus Q(Z) ~ f or } ar g e At criticality the 

correlation length diverges and 

]im R x (U) = l"Q(0) • (36) 

When the system is governed by a disordered fixed point to — 0. When the system is 
governed by a pure fixed point u = {^j . Thus the two possible behaviors for 1 <C £ <C / 



arc 



-d 

, , , for a random fixed point 
Rx(U)~\ , (37) 

1 1 £ v tor a pure fixed point 



i 

e 

In an experiment, since generating many samples is impractical, one studies a single large 
sample with a particular realization of the quenched disorder of size I. For any £ the value of 
X measured in the sample is a sampling from a probability distribution with relative variance 
Rx(£,l)- Thus R x (£,l) controls the deviation of X from the many samples average. If the 
system is governed by a random fixed point, as the correlation length is increased, R x 

increases as ~ (ij . X behaves like the average of (l) independent measurements on 
regions of size £ d . The variance of these measurements does not decrease as £ increases; it 
is constant. On the other hand if the disorder is irrelevant and the system is governed by 
a pure fixed point with a < 0, as the correlation length is increased, R x increases more 

mildly as ~ In this case too, X behaves like the average of (l) measurements 

on regions of size £ d . However, as £ increases, the variance of these measurements decreases 
<-— 

as ~ 4 v ■ 

We have verified that for a disordered system governed by a random fixed point (ST c (l)) 2 
does not scale as ~ l~ d , but rather that (6T c (l)) 2 ~ This is an important result, similar 
to the situation in the purely geometric percolation problem P5| |. Recently Pazmandi, 



Scalettar and Zimanyi |12[ claimed that the bound v > 2/d, which was supposed to hold 
for disordered systems |pTj[ , is not valid. As they show, if in systems violating this bound 
one would have (<5T C (7)) 2 ~ l~ d [our equation (|5])], then simulations at T c °° would not be 
able to capture the true critical exponents [ 12| . In fact in [[Tj|] (^|) is termed "the most likely 



scenario" and the conclusion drawn is that "self averaging breaks down" . However, studies 
of percolation |j45l , our results, and those of AH JTJ| imply the contrary. (<5T C (7)) 2 ~ l~~ [our 



equation (^)], and therefore simulations at T c °° are able to capture the true critical exponents 
even if v < 2/d. This also becomes evident by examining the finite size scaling theory of 
for [Xj(T c °°)], assuming that (^) holds versus the consequences of (|8[). 

We've shown that fluctuations in X^ at T c °° are predominantly due to fluctuations in 
5Ti = T c °° — T c (i,l), and that these fluctuations can be dramatically reduced by measuring 



28 



Xi at T c (i, I). This suggests that using the histogram method to obtain Xi(T c (i, /)) for each 
sample might be a better strategy for Monte Carlo studies than the current strategy of 
studying Xj(T c °°). It was also shown that to a good approximation, fluctuations of Xi close 
to criticality can be accounted for by the finite size scaling form d32|) . We believe that a more 
extensive study of the finite size scaling of sample to sample fluctuations is both feasible 
and desired. 

One of the surprising results of this work is the difference found between the p = 0.6 
model with grand canonical disorder and the c = 0.6 model with canonical disorder. Our 
results indicate that for p = 0.6 and c = 0.6 V#- max scales as l~ 2yt and that V^ max (p = 
0.6)/V/^ max (c = 0.6) ~ 3.26. This is apparently the reason why for these two types of 
disorder Rx tends as I — > oo to different constants. On the other hand we did not find any 
difference in the scaling exponents of the two types of disorder. 
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